
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The impact of divorce laws on the equilibrium in the marriage market.
% Reynoso
% April 2024
%
% This file inputs outputs from the model estimated under MCD
% and the model simulated under UD and constructs marriage market
% and family life cycle outcomes.
%  
% Data: PSID 1968-1992
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------------------ Preliminaries ----------------------------------%

%--- Working directory
cd ([replication_location,'\Outputs']);

%--- Path for estimates
estimates_dir = [replication_location,'\Outputs\MCD\'];

%--- Path to store
output_dir = [replication_location,'\Outputs'];

%--- Load outputs
load([estimates_dir, 'outputs_mcd']);
params=X;
load([output_dir, '\MCD\moments_mc']); 
load([output_dir, '\UD\moments_ud']); 
load([output_dir, '\MCD\U_mc']); 
load([output_dir, '\UD\U_ud']); 
load([output_dir, '\UD\equil_ud']); 

%------------------------ UD full model outcomes -------------------------%

%--- Organize moments
model_moments_bymkt_mc = reshape(model_moments_mc(1:168), size(model_moments_bymkt_ud,1),size(model_moments_bymkt_ud,2));
mf_f_mc = model_moments_bymkt_mc(1:9,:); 
mf_f_ud = model_moments_bymkt_ud(1:9,:); 
sahw_mc = model_moments_bymkt_mc(25:33,:); 
divorce_mc = model_moments_bymkt_mc(34:42,:); 
sahw_ud = model_moments_bymkt_ud(25:33,:); 
divorce_ud = model_moments_bymkt_ud(34:42,:); 
sorting_f_mc = zeros(n_g+1,4);
sorting_f_ud = zeros(n_g+1,4);
for mkt=1:4
sorting_f_mc(:,mkt)=[(mf_f_mc(1,mkt)+mf_f_mc(5,mkt)+mf_f_mc(9,mkt))/sum(mf_f_mc(:,mkt)) ;
    mf_f_mc(:,mkt)./sum(mf_f_mc(:,mkt))];
end
for mkt=1:4
sorting_f_ud(:,mkt)=[(mf_f_ud(1,mkt)+mf_f_ud(5,mkt)+mf_f_ud(9,mkt))/sum(mf_f_ud(:,mkt)) ; 
    mf_f_ud(:,mkt)./sum(mf_f_ud(:,mkt))];
end

%--- Correlation
N_f = 10000; 
NHhs_mc = round(sorting_f_mc(2:10,:)*N_f);
NHhs_ud = round(sorting_f_ud(2:10,:)*N_f);
Tf=[1 1 1 2 2 2 3 3 3]; 
Tm=[1 2 3 1 2 3 1 2 3];
%mkt 1
wcouples_mc1 = repelem(Tf,NHhs_mc(:,1));
mcouples_mc1 = repelem(Tm,NHhs_mc(:,1));
wcouples_ud1 = repelem(Tf,NHhs_ud(:,1));
mcouples_ud1 = repelem(Tm,NHhs_ud(:,1));
%mkt2
wcouples_mc2 = repelem(Tf,NHhs_mc(:,2));
mcouples_mc2 = repelem(Tm,NHhs_mc(:,2));
wcouples_ud2 = repelem(Tf,NHhs_ud(:,2));
mcouples_ud2 = repelem(Tm,NHhs_ud(:,2));
%mkt3
wcouples_mc3 = repelem(Tf,NHhs_mc(:,3));
mcouples_mc3 = repelem(Tm,NHhs_mc(:,3));
wcouples_ud3 = repelem(Tf,NHhs_ud(:,3));
mcouples_ud3 = repelem(Tm,NHhs_ud(:,3));
%mkt4
wcouples_mc4 = repelem(Tf,NHhs_mc(:,4));
mcouples_mc4 = repelem(Tm,NHhs_mc(:,4));
wcouples_ud4 = repelem(Tf,NHhs_ud(:,4));
mcouples_ud4 = repelem(Tm,NHhs_ud(:,4));

corr_mc1 = corrcoef(wcouples_mc1, mcouples_mc1);
corr_ud1 = corrcoef(wcouples_ud1, mcouples_ud1);
corr_mc2 = corrcoef(wcouples_mc2, mcouples_mc2);
corr_ud2 = corrcoef(wcouples_ud2, mcouples_ud2);
corr_mc3 = corrcoef(wcouples_mc3, mcouples_mc3);
corr_ud3 = corrcoef(wcouples_ud3, mcouples_ud3);
corr_mc4 = corrcoef(wcouples_mc4, mcouples_mc4);
corr_ud4 = corrcoef(wcouples_ud4, mcouples_ud4);

change_corr_bymkt = [(corr_ud1(2,1)-corr_mc1(2,1))/corr_mc1(2,1) (corr_ud2(2,1)-corr_mc2(2,1))/corr_mc2(2,1) (corr_ud3(2,1)-corr_mc3(2,1))/corr_mc3(2,1) (corr_ud4(2,1)-corr_mc4(2,1))/corr_mc4(2,1)];
change_corr =[change_corr_bymkt*[0.23 0.23 0.36 0.18]' change_corr_bymkt];
corr_bymkt_mc=[corr_mc1(2,1);corr_mc2(2,1);corr_mc3(2,1);corr_mc4(2,1)];
corr_bymkt_ud=[corr_ud1(2,1);corr_ud2(2,1);corr_ud3(2,1);corr_ud4(2,1)];

%--- Assortative matching measures
fraction_married_mc = sum(mf_f_mc);
hsm=(mf_f_mc(1,:) + mf_f_mc(4,:) +mf_f_mc(7,:))./fraction_married_mc;
scm=(mf_f_mc(2,:) + mf_f_mc(5,:) +mf_f_mc(8,:))./fraction_married_mc;
cpm=(mf_f_mc(3,:) + mf_f_mc(6,:) +mf_f_mc(9,:))./fraction_married_mc;
married_men_byeduc_bymkt_mc = [hsm;scm;cpm];
hsf=(mf_f_mc(1,:) + mf_f_mc(2,:) +mf_f_mc(3,:))./fraction_married_mc;
scf=(mf_f_mc(4,:) + mf_f_mc(5,:) +mf_f_mc(6,:))./fraction_married_mc;
cpf=(mf_f_mc(7,:) + mf_f_mc(8,:) +mf_f_mc(9,:))./fraction_married_mc;
married_women_byeduc_bymkt_mc = [hsf;scf;cpf];

clear hsm scm cpm hsf scf cpf

fraction_married_ud = sum(mf_f_ud);
hsm=(mf_f_ud(1,:) + mf_f_ud(4,:) +mf_f_ud(7,:))./fraction_married_ud;
scm=(mf_f_ud(2,:) + mf_f_ud(5,:) +mf_f_ud(8,:))./fraction_married_ud;
cpm=(mf_f_ud(3,:) + mf_f_ud(6,:) +mf_f_ud(9,:))./fraction_married_ud;
married_men_byeduc_bymkt_ud = [hsm;scm;cpm];
hsf=(mf_f_ud(1,:) + mf_f_ud(2,:) +mf_f_ud(3,:))./fraction_married_ud;
scf=(mf_f_ud(4,:) + mf_f_ud(5,:) +mf_f_ud(6,:))./fraction_married_ud;
cpf=(mf_f_ud(7,:) + mf_f_ud(8,:) +mf_f_ud(9,:))./fraction_married_ud;
married_women_byeduc_bymkt_ud = [hsf;scf;cpf];

%//random matchings

%-MC
rm_mc_mkt1 = married_women_byeduc_bymkt_mc(:,1)*married_men_byeduc_bymkt_mc(:,1)';
rm_mc_mkt2 = married_women_byeduc_bymkt_mc(:,2)*married_men_byeduc_bymkt_mc(:,2)';
rm_mc_mkt3 = married_women_byeduc_bymkt_mc(:,3)*married_men_byeduc_bymkt_mc(:,3)';
rm_mc_mkt4 = married_women_byeduc_bymkt_mc(:,4)*married_men_byeduc_bymkt_mc(:,4)';

%-UD
rm_ud_mkt1 = married_women_byeduc_bymkt_ud(:,1)*married_men_byeduc_bymkt_ud(:,1)';
rm_ud_mkt2 = married_women_byeduc_bymkt_ud(:,2)*married_men_byeduc_bymkt_ud(:,2)';
rm_ud_mkt3 = married_women_byeduc_bymkt_ud(:,3)*married_men_byeduc_bymkt_ud(:,3)';
rm_ud_mkt4 = married_women_byeduc_bymkt_ud(:,4)*married_men_byeduc_bymkt_ud(:,4)';

%//actual matchings

%-MC
am_mc_mkt1 = reshape(sorting_f_mc(2:10,1),3,3);
am_mc_mkt2 = reshape(sorting_f_mc(2:10,2),3,3);
am_mc_mkt3 = reshape(sorting_f_mc(2:10,3),3,3);
am_mc_mkt4 = reshape(sorting_f_mc(2:10,4),3,3);

%-UD
am_ud_mkt1 = reshape(sorting_f_ud(2:10,1),3,3);
am_ud_mkt2 = reshape(sorting_f_ud(2:10,2),3,3);
am_ud_mkt3 = reshape(sorting_f_ud(2:10,3),3,3);
am_ud_mkt4 = reshape(sorting_f_ud(2:10,4),3,3);

%//difference
ggks_mc_mkt1 = am_mc_mkt1 - rm_mc_mkt1;
ggks_mc_mkt2 = am_mc_mkt2 - rm_mc_mkt2;
ggks_mc_mkt3 = am_mc_mkt3 - rm_mc_mkt3;
ggks_mc_mkt4 = am_mc_mkt4 - rm_mc_mkt4;
ggks_ud_mkt1 = am_ud_mkt1 - rm_ud_mkt1;
ggks_ud_mkt2 = am_ud_mkt2 - rm_ud_mkt2;
ggks_ud_mkt3 = am_ud_mkt3 - rm_ud_mkt3;
ggks_ud_mkt4 = am_ud_mkt4 - rm_ud_mkt4;
ggks_mc_bymkt = zeros(3,3,4);
ggks_mc_bymkt(:,:,1)=ggks_mc_mkt1;
ggks_mc_bymkt(:,:,2)=ggks_mc_mkt2;
ggks_mc_bymkt(:,:,3)=ggks_mc_mkt3;
ggks_mc_bymkt(:,:,4)=ggks_mc_mkt4;
ggks_ud_bymkt = zeros(3,3,4);
ggks_ud_bymkt(:,:,1)=ggks_ud_mkt1;
ggks_ud_bymkt(:,:,2)=ggks_ud_mkt2;
ggks_ud_bymkt(:,:,3)=ggks_ud_mkt3;
ggks_ud_bymkt(:,:,4)=ggks_ud_mkt4;

%//Ratio of traces
delta_bymkt_mc=[trace(am_mc_mkt1)/trace(rm_mc_mkt1); trace(am_mc_mkt2)/trace(rm_mc_mkt2); trace(am_mc_mkt3)/trace(rm_mc_mkt3); trace(am_mc_mkt4)/trace(rm_mc_mkt4)];
delta_bymkt_ud=[trace(am_ud_mkt1)/trace(rm_ud_mkt1); trace(am_ud_mkt2)/trace(rm_ud_mkt2); trace(am_ud_mkt3)/trace(rm_ud_mkt3); trace(am_ud_mkt4)/trace(rm_ud_mkt4)];

%--- Gains from marriage
EVM_mc_bymkt = log(U_mc_tot_bymkt) ;
EVMf_mc_bymkt = [EVM_mc_bymkt([1 4 7],1,1) EVM_mc_bymkt([1 4 7],1,2) EVM_mc_bymkt([1 4 7],1,3) EVM_mc_bymkt([1 4 7],1,4)];
EVM_ud_bymkt = log(U_ud_tot_bymkt)  ;
EVMf_ud_bymkt = [EVM_ud_bymkt([1 4 7],1,1) EVM_ud_bymkt([1 4 7],1,2) EVM_ud_bymkt([1 4 7],1,3) EVM_ud_bymkt([1 4 7],1,4)];
GM_mc = [log(U_mc_tot_bymkt(:,1,:)) - U_mc_bymkt(:,2,:)  log(U_mc_tot_bymkt(:,2,:)) - U_mc_bymkt(:,4,:)];
GMf_mc = [GM_mc([1 4 7],1,1) GM_mc([1 4 7],1,2) GM_mc([1 4 7],1,3) GM_mc([1 4 7],1,4)];
GM_ud = [log(U_ud_tot_bymkt(:,1,:)) - U_ud_bymkt(:,2,:)  log(U_ud_tot_bymkt(:,2,:)) - U_ud_bymkt(:,4,:)];
GMf_ud = [GM_ud([1 4 7],1,1) GM_ud([1 4 7],1,2) GM_ud([1 4 7],1,3) GM_ud([1 4 7],1,4)];
%// Consumption equivalence
Pi_c_f = (1- exp((EVMf_ud_bymkt-EVMf_mc_bymkt)/9.15))*100;

%--- Save
save('UD\UD_effects_fullequilibrium', 'corr_bymkt_mc', 'corr_bymkt_ud', 'change_corr', 'GMf_mc', 'GMf_ud', 'Pi_c_f',  'delta_bymkt_mc', 'delta_bymkt_ud', 'ggks_mc_bymkt', 'ggks_ud_bymkt', 'sahw_mc', 'sahw_ud', 'divorce_mc', 'divorce_ud' , 'mf_f_mc', 'mf_f_ud');















